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ABSTRACT 

We present new, high-precision Doppler radial velocity (RV) data sets for the nearby K3V star 
HD 219134. The data include 175 velocities obtained with the HIRES Spectrograph at the Keck I 
Telescope, and 101 velocities obtained with the Levy Spectrograph at the Automated Planet Einder 
Telescope (APE) at Lick Observatory. Our observations reveal six new planetary candidates, with 
orbital periods of P = 3.1, 6.8, 22.8, 46.7, 94.2 and 2247 days, spanning masses of Alsin/ = 3.8, 3.5, 

8.9, 21.3, 10.8 and 108 AI 0 respectively. Our analysis indicates that the outermost signal is unlikely to 
be an artifact induced by stellar activity. In addition, several years of precision photometry with the 
TIO 0.8 m automatic photometric telescope (APT) at Eairborn Observatory demonstrated a lack of 
brightness variability to a limit of ^0.0002 mag, providing strong support for planetary-reflex motion 
as the source of the radial velocity variations. The HD 219134 system, with its bright (V = 5.6) 
primary provides an excellent opportunity to obtain detailed orbital characterization (and potentially 
follow-up observations) of a planetary system that resembles many of the multiple-planet systems 
detected by Kepler, and which are expected to be detected by NASA’s forthcoming TESS Mission 
and by ESA’s forthcoming PLATO Mission. 

Subject headings: Planets and satellites: detection. Planetary systems. Stars: individual (HD 219134), 
Techniques: Radial Velocities 


1. INTRODUCTION 

Extrasolar planets are a source of substantial fasci¬ 
nation and excitement, and with thousands of exam¬ 
ples now known, the statistics of the global distribu¬ 
tion are coming into focus. Goncentrations of planets - 
well-delineated populations in the mass-period diagram - 
have been evident for some time, with super-Earths, hot 
Jupiters and longer-period eccentric giant planets form¬ 
ing groupings that, while distinct, are of effectively still 
unknown province. 

An apt analogy can be drawn with the gradual dis¬ 
covery of the orbital distributions of asteroids within our 
own Solar System. Once minor planets had been discov¬ 
ered in significan t quantity, clear structures such as the 
Kirkwood Gaps (jKirkwood 1866) became progressively 
more apparent, although their origin remained mysteri¬ 
ous. It is sobering to remark that even after more than 
two centuries, and the identification of the gaps as aris¬ 
ing from resonant dynamics, the formation and evolution 
of the asteroids remain topics of active research. 

In this paper, we report that our set of 276 velocities 
for HD 219134 (including 138 measurements with 2-hour 
binning obtained from long-term Keck planet surveys, 37 
measurements with 2-hour binning obtained from spec- 
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tra taken at Keck by the NASA QOl Progranj^ and 101 
measurements with 2-hour binning made with the Auto¬ 
mated Planet Einder (APE) telescope) reveal that this 
star hosts a multi-planet system. 

Indeed, the radial velocities obtained at Keck have, 
since 2010, strongly suggested that HD 219134 is accom¬ 
panied by a multiple-planet system, but the orbital ar¬ 
chitecture at periods P < 100 d was unclear; the observ¬ 
ing cadence at Keck was insufficient to adequately de¬ 
fine the orbital parameters of this rather complex multi¬ 
planet system. We find that the new APE data, however, 
with their high velocity precision and improved observ¬ 
ing cadence, permit much fuller orbital characterizations 
for the planetary candidates. Our best model indicates 
that the star is accompanied by an inner configuration of 
five low-amplitude planets (having radial velocity half¬ 
amplitudes of K=1.9 ms“^, 1.4 ms“^, 2.3 ms“^, 4.4 
ms“^, and 1.8 ms“^, all with orbital periods P < 100 
days). The system also displays a longer-period signal 
with P = 2247 ± 43 d, and Msin(/) = 0.34 ± 0.02 Mj^p, 
which is similar to the mass of Saturn. The presence of 
this outer planet has interesting consequences for current 
planet formation theories. 

Taken as a whole, HD 219134 presents a planetary sys¬ 
tem of substantial scientific interest. Its retinue of multi¬ 
ple super-Earth category planets is highly reminiscent of 
many of the systems discovered by Kepler, albeit in as¬ 
sociation with a star that is thousands of times brighter 
than the median star in the Kepler catalog. 

The plan for this paper is as follows: in §2, we review 
what is currently known about HD 219134. In §3, we 
provide a pro-forma review of our Doppler technique, as 

® NASA program: “TPF Preparatory Science: Low Mass Short- 
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Fig. 1.— HR diagram with HD 219134’s position indicated as 
a small open circle. Absolute magnitudes, M, are estimated from 
V-band apparent magnitudes and Hipparcos distances using M = 
V + 5 log]^Q(d/10 pc). All 956 stars in our catalog of radial velocity 
measurements for which more than 20 Doppler measurements exist 
are shown, color-coded by their B-V values. 

well as an up-to-date report on the performance and re¬ 
cent results obtained with the APF telescope and Levy 
spectrograph. In §4, we discuss our radial velocity ob¬ 
servations for HD 219134, and the 6 planet model that 
we use to interpret the velocity variations exhibited by 
the star. In §5 we discuss our photometric observations 
of the system. In §6 we very briefly assess the system in 
the light of current theories regarding planet formation 
and conclude. 


2. STELLAR PARAMETERS 

HD 219134 (HR 8832; GJ 892; HIP 114622 ) is lo¬ 
cated high in the northern sky (RA = -1-23:13:17, DEC 
+57:10:06). As a bright {V = 5.57), nearby (d = 6.55 pc) 
K-type main sequence dwarf of naked-eye visibility, it has 
long been of interest as a potential planet-bearing star. 
It was among the original 23 UBC Precise Radial Veloc¬ 
ity program star s observed with the CFHT starti ng in 
the early 1980s (Walker et al. 1995 Walker 12012), and 
it was an early target of interest at Keck, i'he first of 
our 138 Keck velocity measurements dates to JD 2450395 
(November, 1996). HD 219134 is currently the 99th near¬ 
est known stellar systenQ and as a consequence, any 
planetary system that it harbors would rank among the 
ten closest known systems (a plurality of which orbit 
much dimmer M-dwarf primaries). Among stars known 
to harbor planets with masses Mpi sin (i) < lOMjup, only 
the Sun, Alpha Cen B, 61 Virginis, and HD 20794 have 
brighter V magnitudes. 

Figure shows HD 219134’s position on a color- 
magnitude diagram containing all of the stars in the 
current Lick-Carnegie Keck database that have accumu¬ 
lated more than 20 Doppler measurements, and empha- 
sizes the star’s entire ly ordinary Main Sequence location. 
Takeda et al. (2007) derive mass and radius estimates 
of MIMq — 0.794 and R/Rc^ = 0.77 , along with an 
age, r = 12.9 Gyr. Valenti & Fischer (2005) measure 
vsm{i) = 1.8kms“^, which, if we assume equator-on ge- 
o metry and a radius Rj R(^ = 0.77, implies Prot ^ 20 d. 

|Tanner et al.| (|2010) used high-contrast imaging of 
HD 219134 with PALO A /PHARO on the Palomar 200- 
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Fig. 2.— The median ^-index values and dispersions of the 
(S-index measurements for the stars in the current Keck sample. 
HD 219134 is shown in red. The size of the points is proportional 
to the number of observations. 


inch telescope to identify three brown-dwarf companion 
candidates at 5 ^7", 5 ^10", and S ^11" separa¬ 
tions. These candidates, however, were determined from 
archiv al second-epoch observations to be background 
stars. Eggleton & Tokovinin (2008), in a multiplicity 
survey of bright nearby stars, report that HD 219134 has 
a V=9.4 companion at S ^106.6" separation. If physi¬ 
cally bound to the primary, this companion would be a 
low-mass M-dwarf with a projected separation d ^ 700 
light years, and an orbital period P ^ 20Kyr. Depend¬ 
ing on the orbital configuration, such a companion might 
be capable of exerting non-negligible long-term gravita¬ 
tional perturbations on the HD 219134 system, and so 
follow-up work to determine the status of the asterism is 
warranted. 

In spite of its potentially great age, HD 219134 does 
show indications of stellar activity. It is listed as a “Flare 
Star” in the SIMBAD Database, and both its median 5'- 
index value and the standard deviation of its individual 
5'-index measurements from our Keck spectra are higher 
than those of t he main locus of stars in o ur Keck survey 
(see Figure 1^. [l saacson fc Fischer| ( |2Q1Q[ ) report a stellar 
velocity jitter of 1.57ms“^ for HD 219134, and this rel¬ 
atively low value is corroborated by the analysis of this 
paper (Section |^. We do find, however, that the ra¬ 
dial velocities for the star are potentially correlated with 
stellar activity over the decade-plus time baseline of our 
observations. A periodogram of the 5'-index values (in¬ 
cluding measurements at all of our Keck epochs, and at 
all APF epochs for which photoelectron counts in the I 2 
region of the spectrum exceed N = 25,000) is shown in 
the top panel of Figure There is a significant peak 
in this periodogram at P ^ 3300 days. This period is 
greater than and district from the P ^ 2300 day period¬ 
icity that is present in the Doppler Velocity data for this 
star. Figure also shows a correlation plot of the ra¬ 
dial velocity observations and their P-index values. The 
peak at zero lag (observation record) indicates that there 
is correlation between the long-period signal in the radial 
velocities and the P-index values, as manifested by the 
long period of decline in both time series. It is there¬ 
fore possible that some of the observed velocity variation 
can be attributed to surface activity. Caution is always 
warranted in interpreting long-term RV variations. 

The Mt. Wilson P-index measures the ratio of flux 
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Fig. 3.— Top Panel shows a periodogram of the Mt. Wil¬ 
son ^-index values associated with our Keck and APF spectra of 
HD 219134. Bottom Panel shows a correlation plot for RV data 
points and their associated (S-index values. The shaded area marks 
the 95% confidence interval for the Pearson correlation coefficient, 
estimated using sets of white noise data. 


from lA bins surrounding the line centers of the Ca II 
H& K lines (at 3968.47A and 3933.66A), as compared to 
two broader 25A bandpasses lying 250A to either side of 
the Ca II H& K line location (Duncan et al. 1991). In 
the standard picture, an increase in the 5'-index, whose 
flux emerges from above the mean photospheric depth of 
the star, corresponds with an increase in spot activity on 
the stellar surface. Spots, in turn, suppress convection 
in their vicinity, which decreases the overall convective 
blueshift of the star, leading to the expectation of a corre¬ 
lation with the Doppler velocity of the star. A star with 
a magnetic cycle that modulates the number of spots can 
therefore present a long-term Doppler trend with an am¬ 
plitude and periodicity that mimics the Keplerian signal 
from a distant planet (Dumusque et al. 2011). 

The upper panel of Figure H charts the velocity mea¬ 
surements taken with the Ke^ and the APF telescopes 
(with the median value for each data set subtracted) 
against the corresponding 5'-index measurements. While 
the strength of the overall positive correlation is indi¬ 
cated by a linear fit to the data, the color-coded time¬ 
ordering of the points indicates that the correlation was 
much stronger during epochs from 2000 through approx¬ 
imately 2010. Indeed, as is indicated by the middle panel 
of Figure the correlation has reversed sign from 2012 
through present, and a weak negative correlation (having 
less than 1 a significance) is present in the recent APF 
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Fig. 4. — Top Panel shows a scatter plot of the radial veloc¬ 
ity data for our Keck and APF observations and their associated 
(S-index value. Each point is marked according to the time of ob¬ 
servation. The best linear fit is shown as a solid line. Middle 
Panel shows the same data, faceted in time to emphasize the time- 
dependent correlation between the radial velocities and the S-index 
activity indicator. Each panel contains the same number of points. 
Bottom Panel shows the same data, faceted by dataset. 


observations. Figure!^ permits comparison of the time 
evolution of the S-ina^ values to the corresponding RV 
observations. 

The unusual features in the time development of the 
RV - 5'-index correlation imply that considerable caution 
must be exercised in interpreting the source of the multi¬ 
year periodicity that is present in the Doppler velocity 
time series. Our candidate planetary signal could, for 
example, be produced by the superposition of the stellar 
magnetic activity cycle and a giant planet on a Keple¬ 
rian orbit. Further monitoring, accompanied by detailed 
analysis, is clearly required. 
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Fig. 5. — Plot of the S-index values compared to the correspond¬ 
ing RV observation. Velocity measurements and S-index values are 
shown here for each I 2 spectrum in our database for HD 219134. 
For purposes of RV modeling, we use 2-hour binning. 


TABLE 1 

Stellar parameters for HD 219134 


Parameter Value 


Reference 


Spectral type 
Mv 
V 

B-V 
Mass {Mq) 
Radius (7^©) 
Luminosity (L©) 
Distance (pc) 
^hk 
Age (Gyr) 
[Fe/H] 
Teff (AT) 
^og{g) (cm s 2) 


K3V 

6.46 
5.57 
0.99 

u. / »^_o.022 
0.77 ± 0.02 
0.31 

6.546 ± 0.012 
0.25 

12.46 
0.08 
4913 
4.51 
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Valenti & Fischer (2005) report a number of additional 
spectroscopically derived properties for HD 219134. It 
appears to be somewhat metal-rich in comparison to the 
Sun, with [M/H] =0.09, and individual abundances that 
include [Na/H] = 0.13, [Si/H] = 0.02, [Ti/H] = 0.02 . 
[Fe/H] = 0.12, and [Ni/H] = 0.09. Ramirez et al. (2013) 
report a high oxygen abundances of [O/Hj = O.zJ We 
note that this high value seems discrepant in light of 
the star’s other abundance m easurements, as well as th e 
value [Fe/H] = 0.04 found by (jAllende Prieto et al.|2004|), 
and so must be treated with caution. One could specu- 
late, however, that there might be a connection between 
the possible high abundance of iron and the apparent 
ease with which th e system formed multip l e planets hav¬ 
i ng P < 100 days ( Robinson et al.|[2006 ). Tanner et al. 


(2009) observed HD 219134 at 160/im using the Spitzer 
Space Telescope’s MIPS spectrometer, and found no ex¬ 
cess emission characteristic of a remnant debris disk, in 
keeping with the large apparent age of the star. 

3. RADIAL VELOCITY OBSERVATIONS 



219134 that form the basis of this paper. In accordance 
with long-established practice, Doppler shifts at both 
telescopes are obtained by imprinting an iodine absorp¬ 
tion spectrum on the collected starlight prior to it s inci¬ 
dence on the spectrograph slit ( Butler et al.||I996 ). The 


forest of added I 2 lines generates a stable wavelength cal 
ibration and permits the measurement of the spectrom¬ 
eter point spread function (PSF). For each spectrum so 
obtained, the 5000 A < A < 6200 A region containing a 
sufficient density of I 2 lines is subdivided into 700 indi¬ 
vidual segments of width 2A, with each segment provid¬ 
ing independent measures of the wavelength, the PSF, 
and the Doppler shift. Our reported overall stellar ve¬ 
locity from a given spectrum is a weighted mean of the 
individual velocity measurements. The uncertainty for 
each velocity is the RMS of the individual segment ve¬ 
locity values about the mean divided by the square root 
of the number of segments. This “internal” uncertainty 
represents primarily errors in the fitting process, which 
are dominated by Poisson statistics. The velocities are 
expressed relative to the solar system barycenter, but 
are not referenced to any absolute fiducial point As a 
consequence, the velocity zero-point offset between the 
measurements at the two telescopes must be treated as 
a free parameter. 

For the data set being considered here, there is an 8- 
year gap between the first Keck velocity measurement 
and the second Keck velocity measurement. The Keck 
HIRES CCD was upgraded during the interval between 
the two observations. In our reduction pipeline, we first 
analyze the entire Keck data set using only the spectral 
chunks that are present in both the pre and post-upgrade 
detector CCDs, thereby obviating the need for any ad¬ 
ditional internal velocity offsets. We then reanalyze the 
post-fix spectra, using all the spectral chunks, thereby 
improving the post-fix precision. 

The Automated Planet Finder (APF) is a 2.4m tele¬ 
scope at Lick Observatory. Coupled with the high- 
resolution Levy echelle spectrometer, it was designed 
to detect planets in the liquid water habitable zone of 
their host stars. It works at a typical spectral resolu¬ 
tion of R ^ 110,000 and delivers a peak overall system 
efficiency (fraction of photons striking the tel escope pri¬ 


mary t hat are detected by the CCD) of 15% (Vogt et al 
2014b). Currently, 80% of the telescope’s time is dedi¬ 
cated to the detection of extrasolar planets, and the APF 
has scheduling software capable of making decisions on 
what target to observe based on the ambient atmospheric 
transparency, atmospheric seeing, and moon phase. This 
allows the telescope to operate efficiently throughout the 
year without the need for human supervision. 

Since it began acquiring scientific data in Q2 2013, the 
telescope has contributed to three planetary system dis- 


coveries, all with radial velocitv measurements having a 

1-3 ms ^ level of precision ( 

Vogt et al.|2014a Burt et al. 

2014 

Fulton et al. 

2015 

). ' 

I'he detection of these plan- 

ets, a] 

long with the cand 

[idates described in the current 


paper, indicates that the APF telescope is well-suited 
to the discovery of low-mass planets orbiting low-mass 
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Internal Uncertainty 


Fig. 6. — Histogram showing internal precisions obtained to date 
with APF while observing bright, nearby Main Sequence Stars. 
The subset of measurements taken of HD 219134 are shown in 
green. Observations are unbinned. 



Time [JD] 

RV [ms-i] 

Uncertainty [ms 

Dataset 

1 

2450395.74 

-4.50 

0.50 

KECK 

2 

2453239.05 

-2.14 

0.74 

KECK 

3 

2453301.77 

-13.31 

1.37 

KECK 

4 

2453338.71 

-5.86 

0.53 

KECK 

5 

2453547.10 

1.56 

0.45 

KECK 

6 

2453548.09 

1.39 

0.48 

KECK 

7 

2453549.12 

-0.66 

0.55 

KECK 

8 

2453550.09 

2.79 

0.43 

KECK 

9 

2453551.09 

0.81 

0.38 

KECK 

10 

2453552.05 

-1.88 

0.44 

KECK 

11 

2453571.06 

0.00 

0.48 

KECK 

12 

2453692.77 

6.20 

0.71 

KECK 

13 

2453693.70 

1.81 

0.60 

KECK 

14 

2453693.92 

1.54 

0.57 

KECK 

15 

2453694.70 

1.48 

0.58 

KECK 

16 

2453694.85 

1.99 

0.61 

KECK 

17 

2453695.85 

7.63 

0.59 

KECK 

18 

2453696.71 

6.48 

0.55 

KECK 

19 

2453713.68 

-1.53 

0.64 

QOl 

20 

2453713.68 

-1.19 

0.65 

QOl 


stars. 

The APF has consistently achieved internal velocity 
precision of order a < 2nis“^ on bright (e.g. V ^ 7) 
stars. Figure 1^ is a histogram showing the internal pre¬ 
cisions obtained to date for scientific target stars in our 
Doppler velocity program. The median internal (photon 
shot noise) precision is 1.18 ms“^, and 2485 of the mea- 
s urements (binned at 2-hour cadence) have a < 2ms“^ 


TABLE 2 

Radial Velocity observations (sample) 


Fulton et al. (2015) report use of radial velocity data 
obtained using the APF telescope to characterize a 
multiple-planet system orbiting HD 7924. In the course 
of their analysis, they investigated correlations with en¬ 
vironmental parameters such as air pressure and CCD 
temperature, and find improvements up to a factor of 
two in the RMS of APF velocities when the data is de¬ 
trended against these factors. 

In light of these results, we initiated a similar set of 
experiments. We compared the velocity values for the 
RV standard stars HD 185144, HD 9407 and HD 10700 
against all of the environmental parameters stored for 
each observation. We detect no notable correlation be¬ 
tween the Doppler velocity measurements and any of 
these factors. When fitting linear trends to the veloci¬ 
ties as a function of air pressure, we see the RMS change 
from 1.84 ms“^ to 1.81 ms“^ for HD185144, 2.13 ms“^ 
to 2.11 ms“^ for HD 10700 and 2.46 ms“^ to 2.33 ms“^ 
for HD 9407. Similarly, when fitting against the CCD 
temperature we see changes in the RMS of 1.84 ms“^ 
to 1.71 ms“^, 2.13 ms“^ to 2.13 ms“^ and 2.46 ms“^ 
to 2.42 ms-i for HD 185144, HD 10700 and HD 9407, 
respectively. 

In our view, it is likely that the difference in correlation 
strengths stems from the different instrument focusing 
procedures and data reduction pipelines used in the sep¬ 
arate analyses of the HD 7924 and HD 219134 spectra. 
Due to the lack of evident trends, we have elected not to 
decorrelate our data set against environmental variables. 

Table presents the complete set of our RV obser¬ 
vations for HD 219134. The RV coverage spans ap¬ 
proximately 19 years of monitoring over 276 (two-hour 
binned) measurements. The median internal uncertainty 
for our observations is ~ 0.75 ms“^, and the peak- 
to-peak velocity is ^ 31.3 ms“^. The velocity scatter 
around the average RV is ^ 5.7 ms“^. Observations 



Fig. 7.— Radial velocity observations for HD 219134. 

marked as “Keck” are HIRES velocities from spectra 
obtained by the Lick-Carnegie Exoplanet team, or, in 
some cases, from publicly available archi ved spectra ob¬ 
taine d by the California Planet Survey (Howard et al. 
2010). Observations marked “QOl” are velocities com- 
puted using our pipeline from archival HIRES spectra 
from the 2005 NASA program: “TPE Preparatory Sci¬ 
ence: Low Mass Short-Period Companions to TPE Tar¬ 
get Stars”, to Principal Investigator W. Cochran. Obser¬ 
vations marked “APE” are from spectra obtained with 
the APE telescope and Levy Spectrometer. 

4. KEPLERIAN SOLUTION 

Eigurej^ shows the RV measurements after binning to 
two-hour increments. We again note that there is a sin¬ 
gle Doppler measurement at BJD 2450395.74 (Novem¬ 
ber 8th, 1996) followed by a nearly eight-year gap to the 
next measurement at BJD 2453239.05 (August 21, 2004). 
The single early epoch point is useful for cementing the 
lack of any apparent long-term large-scale Doppler ve¬ 
locity trend. The top panel of Eigurej^ shows the error- 


weighted, normalized Lo mb-Scargle periodogram (Zech- 
meister & Kiirster 2009). The three horizontal lines m 
the plot represent ditterent levels of false alarm proba- 
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Fig. 8. — Top panel: Error-weighted Lomb-Scargle periodogram 
for HD 219134. False-alarm probability levels are shown at the 
10%, 1% and 0.1% level. Bottom panel: Spectral window function. 


bility (FAP; 10%, 1% and 0.1%, from bottom to top, 
respectively). The FAPs were computed by scrambling 
the dataset 100,000 times and sampling the periodogram 
at 100,000 frequencies, in order to determine the proba¬ 
bility that the power at each frequ ency could be exceeded 
by chance (e.g. Marcy et al.||2005). The bottom panel of 
Figure^ shows the spectral window, displaying the usual 
peaks due to observational cadence, arising fro m the side¬ 
real and solar d ays, and from the solar year (Dawson & 


Fabry cky 


2010 ). 


We tit the radial velocities with a Keplerian model with 
a vector of parameters consisting of the orbital el¬ 
ements (period, mass, mean anomaly, eccentricity and 
longitude of pericenter for each planet) and vertical off¬ 
sets for each dataset (to account for differences in the 
zero point among datasets). Each radial velocity mea¬ 
surement Vi^ taken at time is represented as 


Vi — V(ti^ 0^ 6i s 


3 ’ 


( 1 ) 


where V{ti^ 0) is the predicted velocity, and is normally 
distributed with variance e? (fixed, and corresponding to 
the formal uncertainties quoted by the observer). The 
term Sj accounts for additional sources of scatter (e.g. 
underestimated measurement errors, stellar jitter, and 
other astrophysical sources of RV variation), modeling 
the residual noise in each 7 -t h dataset as n ormally dis¬ 


tributed with variance (e.g. Gregory 2005). Therefore, 


scatter from the model is modeled with a Gaussian dis¬ 
tribution of variance e? + 5^. The best-fit parameters are 


derived by optimizing the log-likelihood of the model: 

No 


logC=-l 


+ E log(27r) 


i=l 


where 


No 


= E Ni - Vif! (e? + 4) ■ 


, ( 2 ) 


(3) 




In order to derive the starting values of the parameters, 
we fit our data by removing peaks in the Lomb-Scargle 
periodogram of the residuals. For each W-planet model, 
we compute bootstrapped periodograms and investigate 
the strongest peaks with FAP < 10“^. Figure shows 
the periodogram of the residuals for each stage m the fit 
construction. 

A simple Markov-Ghain Monte Garlo a lgorithm 
(^MGMG; e.g. Ford 2005 2006 Gregory 2011), in con¬ 
junction with Equations 1-3 and hat priors on log P, 
log and the other orbital parameters, was used to 
characterize the distribution of the parameters of the 
model, using the best fit the starting point. Eor the 
noise parameters, sj, the corresponding prior is a modi¬ 
fied Jeffrey function p(5j) = [(sj+so) ln(l + Smax/«5o)]; for 
Sj <C So (which we take equal to 0.3ms“^), the function 
is a uniform prior (which includes 0), w hile for s . > g o 
the function is a regular Jeffrey prior ([Gregory 2011). 
The MGMG routine is run until sufficient convergence 
is achieved. The uncertainties in each parameter are re¬ 
ported in Table pT] within square brackets. The marginal 
distribution of each parameter, based on 5 x 10^ samples 
from the MGMG routine, is shown in Eigur e[TT| 

The final best-fit model is shown in Eigure [ToT w ith the 
corresponding orbital elements listed in Table |10] where 
we report median and mean absolute deviationior each 
parameter. The inner 5 planets are fixed on circular or¬ 
bits, since the best-fit Keplerian model with eccentric 
orbits has crossing orbits and is therefore unstable. We 
also note that the gravitational (non-Keplerian) interac¬ 
tion between the inner 5 planets is small, but not com¬ 
pletely negligible, since the planet pair b-c and the planet 
triple d-e-f lie close to mean-motion resonances (2:1 and 
4:2:1 respectively). Eigure 12 shows an orbital diagram 
of the system. 

As a final test to assess the quality of the 6-planet 
model, we use a cross-validation algorithm on the data, 
comparing an Wplanet model with an {N — 1)-planet 
model. In the “leave-one-out” flavor used for the present 
paper, we divide the full dataset of observations into 
a training set of A'o ~ 1 observations and a testing set 
of a single observation, rotated among all observations; 
each training set is used to derive a new fit. The like¬ 
lihood of the prediction made combined from each test¬ 
ing set (log>Ccv; higher is better) measures the predictive 
power of the model, and is sensitive to both underfitting 
and overfitting; for instance, a lower (worse) log>Ccv for 
a model with a larger set of parameters is indicative of 
overfitting. Table reports the values of log/icv Each 
fit is strictly better (higher likelihood) than the previ¬ 
ous, suggesting (but not conclusively proving) that the 
6-planet Keplerian model is not overfitting the data. 

Eigure [^illustrates that the distribution of the resid¬ 
uals is very nearly normal, which suggests that we are 
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[3] Residuals for 2-planet fit (g, e) 
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Period [d] 


[4] Residuals for 3-planet fit (g, e, d) 


10 100 1000 10000 
Period [d] 



[6] Residuals for 5-planet fit (g, e, d, b, f) 



100 1000 
Period [d] 


[7] Residuals for 6-planet fit (g, e, d, b, f, c) 



10 100 1000 
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Fig. 9.— Panels (l)-(7): Lomb-Scargle periodograms computed for each model (0-planet, 1-planet, 2-planet, 3-planet, 4-planet, 5-planet 
and 6-planet fits). 




























































































HD 219134b 

HD 219134c 

HD 219134d 

HD 219134e 

HD 2191341 

HD 219134g 

P [days] 

3.0931 [0.0001] 

6.7635 [0.0006] 

22.805 [0.005] 

46.71 [0.01] 

94.2 [0.2] 

2247 [43] 

A4sin(i) [Aljup] 

0.012 [0.001] 

0.011 [0.002] 

0.028 [0.003] 

0.067 [0.004] 

0.034 [0.004] 

0.34 [0.02] 

M [deg] 

57 [20] 

78 [27] 

263 [20] 

277 [11] 

107 [35] 

209 [56] 

e 

0 

0 

0 

0 

0 

0.06 [0.04] 

w [deg] 

0 

0 

0 

0 

0 

215 [50] 

K [ms-l] 

1.9 [0.2] 

1.4 [0.2] 

2.3 [0.2] 

4.4 [0.2] 

1.8 [0.2] 

6.1 [0.3] 

a [AU] 

0.0384740 [8 x 10“^] 

0.064816 [4 X 10-®] 

0.14574 [2 X 10-®] 

0.23508 [4 X 10-®] 

0.3753 [0.0004] 

3.11 [0.04] 

^peri [JD] 

2449999.5 [0.2] 

2449998.5 [0.5] 

2449983 [1] 

2449964 [1] 

2449972 [9] 

2448725 [356] 


QOl 5noise [^^ ^] 
KECK Snoise [ms-^i 
APE Snoise [mS ] 
QOl Av [ms-^] 
KECK Av [ms-i] 
APE Av [ms-^] 

1.1 [0.2] 

2.5 [0.2] 

1.8 [0.2] 

-0.9 [0.6] 

-0.8 [0.2] 

-2.3 [0.6] 

Mi, [Al©] 

0.794 


280.407 

-log/: 

593.311 

RMS [ms-i] 

2.223 

Jitter [mS-1] 

2.038 

Epoch [JD] 

2450000 

Data points 

276 


TABLE 3 

Best-fit 6-Keplerian Model for HD 219134 



RV, Planet e [m/s] RV, Planet d [m/s] RV, Planet c [m/s] RV, Planet b [m/s] 
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Fig. 10.— Best-fit Keplerian model. 











10 



3.0926 3.0930 3.0934 

Period of planet b 



6.761 6.763 6.765 6.767 

Period of planet c 



22.78 22.80 22.82 

Period of planet d 



0.006 0.010 0.014 

Mass of planet b 



0.005 0.010 0.015 0.020 

Mass of planet c 



0.015 0.025 0.035 

Mass of planet d 



Mean anomaly of planet b 



° 0 50 100 150 200 

Mean anomaly of planet c 



200 250 300 350 

Mean anomaly of planet d 



0.06480 0.06482 

Semi-major axis of planet c 



0.14565 0.14575 0.14585 

Semi-major axis of planet d 



1.0 1.5 2.0 2.5 

Semiamplitude of planet b 



0.5 1.0 1.5 2.0 

Semiamplitude of planet c 



1.5 2.0 2.5 3.0 

Semiamplitude of planet d 



46.66 46.70 46.74 

Period of planet e 



0.050 0.060 0.070 0.080 


Mass of planet e 



93.6 94.0 94.4 94.8 


Period of planet f 



0.02 0.03 0.04 0.05 

Mass of planet f 



‘^050 2150 2250 2350 

Period of planet g 



0.25 0.30 0.35 0.40 

Mass of planet g 



240 260 280 300 320 

Mean anomaly of planet e 



° 0 50 100 150 200 

Mean anomaly of planet f 



° 0 50 150 250 350 

Mean anomaly of planet g 



0.2349 0.2351 

Semi-major axis of planet e 



0.3735 0.3745 0.3755 0.3765 

Semi-major axis of planet f 




3.5 4.0 4.5 5.0 5.5 

Semiamplitude of planet e 



1.0 1.5 2.0 2.5 

Semiamplitude of planet f 



4.5 5.0 5.5 6.0 6.5 7.0 7.5 
Semiamplitude of planet g 



0.00 0.10 0.20 0.30 

Eccentricity of planet g 



° 0 50 150 250 350 

Long, of periastron of planet g 



0.5 1.0 1.5 2.0 

Noise, Q01 



2.0 2.5 3.0 3.5 

Noise, KECK 



1.5 2.0 2.5 

Noise, APE 


Fig. 11 


Marginal distributions of the orbital elements, drawn from 5 X 10^ Markov-Chain Monte Carlo samples, 
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Fig. 12. — Left panel: Orbital plot of the 6-planet model. Right panel: Orbital plot of the inner 5 planets. The radius of each point is 
proportional to the square root of its minimum mass. 

log i2cv 


No planets 

-804.34 

1 planet (g) 

-754.59 

2 planets (g, e) 

-699.17 

3 planets (g, e, d) 

-672.69 

4 planets (g, e, d, b) 

-640.45 

5 planets (g, e, d, b, f) 

-614.71 

6 planets (g, e, d, b, f, c) 

-594.11 


TABLE 4 

Cross-validation results 


not overfitting the data (assuming the real distribution 
of the residuals is itself normal). Finally, we note that 
the orbital elements (for signals with P < 100 d) derived 
with the APF data alone are comparable to the orbital 
elements of the fit derived with the combined datasets, 
aside for the 94-day candidate planet, which is detected 
at the lowest significance among the six signals. 

4.1. Risk Assessment 

Figure pA] is a mass-period diagram for planets whose 
detection^ave been publicly announced but whose sta¬ 
tus is currently listed as “controversial” in the exo¬ 
planet.eu database. Nearly 200 objects are plotted; the 
concentration with the region of small K (that is, low 
mass), and P < 100 days is both noticeable and sober¬ 
ing. The use of precision Doppler velocity measurements 
to detect the class of planetary systems that dominate 
the Kepler census is fraught with potential pitfalls. 

The time series of radial velocity measurements for HD 
219134 is no exception, and there are several specific con¬ 
cerns. (1) With half-amplitudes K < 2 ms“^ for three of 
the inner planets (b,c, and d), the signals are weak. (2) 
For planets with low amplitudes, the presence of aliases 
can plague correct interpretation of the periodicities in 
the data. (3) The proximity to mean motion resonances 
for b-c and for d-e-f leads to inconsistencies between N- 
body and Keplerian models of the system. (4) The ec¬ 
centricities of all the inner candidates must be very low 
in order to ensure orbital stability over the long term. 

The residuals periodogram for a 5-Keplerian solution 



Theoretical Quantiles 



Nornnalized residuals 

Fig. 13. — Top panel: Quantile-quantile plot of the residuals 
from the 6-planets model. Perfectly normally distributed residuals 
would fall on the solid line. Bottom panel: Kernel density for the 
normalized residuals (including noise) for each dataset, compared 
to a Gaussian density distribution (dotted line). 

containing the 3.1, 22.8, 46.7, 94, and 2247 day candidate 
periods reveals remaining peaks at P ~ 6.76 days and 
P ^ 1.17 days (Figurefl^. In our analysis, we have so far 
assumed that the 1.1 p^y period is an alias of the 6.76 
day signal. Concerns regarding aliases are heightened for 
the 6.76 day signal as a consequence of the K = 1.4 ± 
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Fig. 14.— Planets listed in the exoplanet.eu catalog as “contro¬ 
versial” . 


30 
25 
20 

0 

I 15 

Q_ 

10 

5 
0 

1 10 
Period [day] 

Fig. 15.— Periodogram for the HD 219134 residuals after sub¬ 
traction of a five-planet solution. There are two power-excess peaks 
around P ~ 1.17 days and P ~ 6.76 days (indicated by the red 
points), one of which is likely the daily alias of the other. 
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Fig. 16.— Periodogram of HD 219134 for the 6.76-day planet 
c (best five-planet solution have been subtracted). The black line 
represents the periodogram of real data sets. The red line repre¬ 
sents the periodogram of an artificial 6.76-d signal + bootstrapped 
residuals to six-planet solution. With artificial signals injected at 
P = 6.76 days, the 1.17-d peak in the observed data sets is fully 
recovered. 


0.2 m RV half-amplitude. To improve confidence that 
the 6.76 day signal represents the true physical period, 
we proceed with the following steps. 

We perform bootstr ap simulations, outlined by |Daw-| 


son & Fabrycky (2010) and adopted in previous analyses 
of aliases in systems containin g multiple low-amplitude 
planets (e.g. Lovis et al. 2011). We generate 500 boot¬ 
strapped Doppler time series by shuffling the residuals 



6.70 6.75 6.80 6.85 6.90 

Period [day] 

Fig. 17.— Periodogram of HD 219134 for the 6.76-day planet 
c (best five-planet solution have been subtracted). The black line 
represents the periodogram of real data sets. The red line repre¬ 
sents the periodogram of an artificial 1.17-d signal + bootstrapped 
residuals to six-planet solution. When injecting an artificial 1.17-d 
signal, the 6.76-d peak in the observed data sets cannot be correctly 
reproduced. 

from the six-planet solution (with either Pq = 6.76 d 
or Pq = 1.17d). For each trial, we inject an artificial 
noiseless sinusoid at the corresponding period, Pq^ and 
compare the resulting periodogram to the observed one. 
Figure and Figure show the results of this exer¬ 
cise, and indicate that the observed periodogram of the 
residuals to the five-planet solution is much more likely 
to be correctly reproduced by a Pe = 6.76 d signal than 
a Pe = 1.17 d signal. 

Additionally, we fit six-Keplerian models to the radial 
velocities using sixth-planet periodicities of Pq = 6.76 d 
and Pq = 1.17d. The 6.76-d signal provides a better fit 
with a reduced = 9-53, and in this case, the 1.17- 
d peak does not remain in the residuals periodogram. 
The six-planet solution with P = 1.17 d yields a worse 
reduced = 10.96. Moreover, the 6.76-d signal is still 
clear in its residuals periodogram. We therefore conclude 
that the peak at P = 6.76 d is a true signal, whereas the 
peak at P = 1.17 d is its one-day alias. 

Another point of concern is the possibility that one (or 
more) of the candidate planet signals are an artifact of 
the stellar rotation period. In §2, we have noted that HD 
219134’s radius and vsm{i) suggest a rotation period of 
order P ^ 20 days, which is close to the period of the 22.8 
day planet candidate, raising the possibility that planet 
“d” is a rotationally modulated spot signal. To check 
whether this might be true, we have looked at differ¬ 
ent epochs of data to determine that the amplitude and 
the period of the signal is constant (a spot signal would 
likely attenuate over a few rotation periods, but could 
then reappear, and would be expected to show significant 
period and frequency drift). We have taken the three 
datasets (Keck, QOl, and APF) and (1) applied the ve¬ 
locity offsets from our best fit, (2) quadrature-augmented 
the noise by the fitted values of as = 1.10ms“^ m/s for 
QOl, aK = 2.50 ms“^ for Keck, and a a = 1.80 ms“^ for 
APF, and (3) sorted the joint data by ascending time 
stamp. We then divided the joint time-sorted data set 
into three segments, each with equal total signal-to-noise 
in its component points. The first data segment thus con¬ 
tains QOl and Keck data, the second segment contains 
Keck and APF data, and the third contains only APF 
data. We subtract out the radial velocity signals from 
the 3.1, 6.7, 46.7, 94.2, and 2247 day candidate planets 
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TABLE 5 

Photometric Amplitudes at HD 219134 Candidate Periods 


Period 

Semi-Amplitude 

RMS Error 

3.093 

0.00025T0.00016 

0.00204 

6.763 

0.00010T0.00016 

0.00205 

22.81 

0.00040T0.00016 

0.00201 

46.72 

0.00025T0.00017 

0.00204 

94.19 

0.00014T0.00017 

0.00202 



Fig. 18.— Differential photometric measurements of HD 2191 34 
using the TIO 0.8m APT at Fairborn Observatory ( |Henry|1999| >. 

using the parameters listed in Table 3. We then compute 
the residuals periodogram for each of the three data seg¬ 
ments. In each case, the P =22.8-day periodicity is the 
location of the strongest peak in the power spectrum of 
the velocity residuals (between P = 10 days and P = 30 
days). This suggests that the P = 22.8-day signal is not 
the product of spot modulation, and that it has been 
present and stable throughout the full time span of our 
observations. 

5. PHOTOMETRIC OBSERVATIONS 


High-precision long-baseline photometric data have 
been acquired for HD 21 9134 with th e TIO 0.8m APT 
at Fairborn Observatory (Henry 1999) in the Stromgren 
b Sz y pass bands. A total of 31o observations were ob¬ 
tained from the 2010 through 2014 observing seasons. 
The two-color observations have been combined to pro¬ 
duce a A(6-|-^)/2 joint-filter time series, which improves 
measurement precision. The time-series is obtained us¬ 
ing the standard quartet observing sequence that com- 
pares the targ et star with a set of three comparison stars 
(Henry||1999). For the observations reported here, the 
three comparison stars (denoted a, b, & c) are: a - HD 
223421, (V = 6.36, B-V=0.408, F2 IV), b - HD 217071 
(V=7.45, B-V=0.368, FI HI), c - HD 215588 (V=6.45, 
B-V=0.430, F5 V), whereas the target star HD 219134 
(5.57, B-V=1.000, K3 V) is denoted star d. 

In reducing the data, all six permutations of differential 
m agnitudes of t he four stars are evaluated. As described 


m 


Henry (1999), only data that survive a cloud filter are 


retained, and the photometry is normalized so that all 
observing seasons have the same mean as the first season. 
We note, however, that there is little, if any, observable 
change in the mean magnitude from year to year. A 3-cr 
filter was applied, which removed six outlying photomet¬ 
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Fig. 19.— Differential photometric measurements of HD 2191 34 
using the TIO 0.8m APT at Fairborn Observatory ( |Henry|1999| >. 


ric points. In keeping with standard procedure, however, 
the outliers were not removed until after the data had 
been phased to each of the planetary periods to be sure 
they were not transit points. 

The full photometric time series is shown in Figure 18 
and a year-by-year breakdown is shown in Figure The 
observations within each year are constant from night to 
night with a mean standard deviation of 0.0014 mag. 
This is the approximate limit of precision for a single 
nightly observation with the TIO APT. The yearly means 
are also constant to a limit of only 0.00055 mag. No sig¬ 
nificant photometric period is found in any year or across 
the entire data set. In addition, we have normalized the 
data so all yearly means are identical and fit least-squares 
sine curves on the candidate planetary periods with the 
results shown in Table 5. None of the candidate plan¬ 
etary periods exhibit significant photometric periodici¬ 
ties. Phase plots on the planetary periods show no sign 
of transits. 

The five inner planets of the HD 219134 system have 
a-priori geometric probabilities of transit of 9.2%, 5.4%, 
2.4%, 1.5% and 0.93% for the P =3.093, 6.763, 22.81, 
46.72, and 94.19 da y periods, respectivel y . As reviewed 
by authors such as Wolfgang & Lopez (2015), super- 
Earth mass planets that are members of systems con¬ 
taining multiple planets with P < 100 d display a very 
large range in radii at given mass, and in expectation, 
often contain ^ 1% of the total planetary mass in hydro¬ 
gen and helium. The presence of these light gasses, in 
turn, generates a substantial contribution to the overall 
planetary radius. If we use the sol ar system mass-radiu s 
relation, Mp/M^ = (Lissauer et al. 2011), 

we expect that the transit depths tor HD 219134 b and 
c will be of order St < 0.1%, which is too small a signal 
for phase-folded long-term ground-based photometry of 
the type reported here, but is within reach if platforms 
such as MOST, Warm Spitzer, or JWST are employed 
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for a targeted transit check. Indeed, systems such as 
HD 219134 form a strong basis for the scientific case of 
the forthcoming CHEOPS Mission, scheduled for launch 
in 2017. 


port from NSF grants AST-0307493 and AST-0908870. 
RPB gratefully acknowledges support from NASA OSS 
Grant NNX07AR40G, the NASA Keck PI program, and 
from the Carneede Institution of Washington. SM ac- 


6. DISCUSSION 

In comparison to our own Solar System, HD 219134 
has an exotic architecture, with at least five super- 
Earth mass planets orbiting with periods of less than 
100 days. Discoveries in recent years, however, have 
i ndicated that suc h systems are surprisingly common 
(Mayor et al.|2009). This result has received strong, and 
i ndeed dramatic con firmation from the Kepler Mission 
(Batalha et al.||2013|), which revealed hundreds of candi- 
date multiple-transiting multiple-planet systems that (at 
least in broad-brush strokes) call to mind HD 219134 b-f. 
This resemblance is underscored by Figure in which 
the HD 219134 planets are shown in conjunction with 
Kepler’s transiting planet candidates. 

A question of substantial interest is whether pla netary 
systems such as HD 219134 are assembled in situ (|Mont- 
gomery fc Langhlin1|2Q091 [Hansen fc: Mnrra 3 ^| 2012 [ |CTm 
ang T Laughlin||2013|), or whether the planets form at 
large distances and then migrate inward to their final 
locations. At present, it is not fully clear how to real¬ 
istically distinguish between the tw o scenarios. Recent 
work by Batygin & Laughlin (2015) has emphasized the 
role of outer giant planets in triggering collisional cas¬ 
cades among planetesimals that can potentially destroy 
systems of super-Earths such as those described in this 
paper. It will therefore be of substantial interest to un¬ 
derstand the nature of the long-term (P > 2000 day) 
periodicity in the radial velocity time series. 

Among the thousands of planetary systems that are 
now known, HD 219134 stands out. The bright primary 
star has demonstrated excellent radial velocity stability 
over two decades of measurement, and, given more data, 
there is a tantalizing possibility of finding additional low- 
mass planets in the system. With the parent star lumi¬ 
nosity estimated at 0.31 Lq, a planet orbiting HD 219134 
at 0.56 AU would receive the same energy flux that the 
Earth receives from the Sun. Such a planet would have 
an orbital period of 167 days. If it had a mass equal to 
that of Earth, its radial velocity half amplitude would 
he K = 14cms“^. Such a signal would be challeng¬ 
ing, but given current projections for the Doppler veloc¬ 
ity technique, almost certainly not impossible to detect. 
Going forward, HD 219134 looks to be an ideal target 
for platforms such as APE, HARPS-N, the APTs, and 
other high-precision Doppler and photometric facilities 
with access to the far northern sky. 
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